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Abstract 

Tsunami-induced coastal floods are relevantly represented 
among the deadliest and costliest disasters ever. In this 
context, several numerical techniques have been studied for 
prevention, protection and planning activities, for both 
catastrophic and more frequent events. An effective solution 
is provided by the SPH method, a mesh-less CFD 
(Computational Fluid Dynamics) technique, particularly 
suitable for flood event modelling. 

This study describes the validation of a 3D "Weakly- 
Compressible" SPH (Smoothed Particle Hydrodynamics) 
model, considering a 3D tsunami test case (Benchmark 1- 
2009 of ISEC-Inundation Science&Engineering Cooperative-). 

The model satisfactorily reproduces the sea level and 
velocity time series, recorded at the monitoring points. The 
results show the reliability and the further potentials of the 
model in representing a 3D generation and propagation of a 
solitary wave with the consequent coastal flood. 

Finally this study seems to represent one of the first 
validations of a 3D SPH model on a tsunami test case. 

Keywords 

SPH; Tsunami; Surface Waves; Floods; Wave Motion; Solitary 
Wave 

Introduction 

Flood events have caused huge damages to human 
health, environment, anthropic structures and 
infrastructures. In particular, tsunami-induced coastal 
floods are heavily represented in both deadliest flood 
events and costliest disasters ever. 

The five deadliest floods, causing over 2.0*10 5 
casualties, are cited in the following. After some years 
of drought, the Central China floods of 1931 (due to 
several months of huge precipitations) caused more 
than 2.5*10 6 dead. Yellow River (China) was interested 
by two huge floods, responsible for more than 9.0*10 5 
casualties in 1887 and more than 5*10 5 dead in 1938. 
Around 2.3*10 5 individuals died because of a dam 
break event (Banqiao, typhoon Nina, China 1975). 


Finally, an Indian Ocean tsunami caused around 
2.3*10 5 casualties in Indonesia (2004). 

On the other hand, the costliest floods ever are 
represented by the following events (damage over 5 
billions USA dollars -USD- each): 

300*10 9 USD: Tohoku earthquake and tsunami 
(2011, Japan), the costliest disaster ever, the 
most damage were due to the tsunami coastal 
flood; 

46*10 9 USD: Thailand riverine and areal floods 
due to monsoons (2011); 

26*10 9 USD. 1998 Yellow river floods; 

7.0*10 9 USD. Armero tragedy (Colombia, 1985, 
pyroclastic flow floods); 

5.0*10 9 USD: 2013 Alberta (riverine and areal) 
floods (Canada). 

In this frame, several numerical techniques have been 
studied in order to simulate flood events (e.g. Le 
Veque et al. 2007), for prevention, protection and 
planning activities. Among the others, SPH (Smoothed 
Particle Hydrodynamics) method is particularly suitable 
for modelling these phenomena, such as coastal floods 
(caused by solitary waves/tsunamis, etc.), dam breaks 
(e.g. Grenier et al. 2009, Di Monaco et al. 2011, 
Amicarelli et al. 2013), riverine and areal floods, urban 
flooding, tank leakages (e.g. Maurel et al. 2009), floods 
caused by landslides or volcanic activities, etc. 

SPH is a mesh-less Computational Fluid Dynamics 
(CFD) technique, whose computational nodes are 
represented by numerical fluid particles. In the 
continuum, the functions and derivatives appearing in 
the governing balance equations of fluid mechanics 
are approximated using convolution integrals over 
smoothing functions. Then for a generic function f the 
SPH integral approximation (<>i) is provided by: 

</>, = I fW<b? (1) 
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where W is a weighting function (called kernel), xo is 
the position of a generic computational point and Vh 
the integration volume called kernel support (a sphere 
of radius 2h in the inner domain, far from boundaries), 
whose characteristic length is h. A generic first order 
derivative along Xi direction can be still computed as 
in (1). After integration by parts, it takes the form: 



J fWn t dx 2 


, , dW 3 

f f dx 

v, dx, 


(2) 


Here a surface integration is considered all over the 
surface (Ah) of the kernel support. This term is 
peculiar and necessary at boundaries, where the 
kernel support is truncated by the fluid domain 
frontiers and SPH truncation errors noticeably grow 
(Amicarelli et al., 2011a and 2011b). The modelling of 
this term can relevantly differentiate SPH models and 
is still object of recent studies (Adami et al. 2012, 
Hashemi et al. 2012, Macia et al. 2012, Mayrhofer et al. 
2013, Ferrand et al. 2013, Amicarelli et al. 2013,). 


SPH technique has several advantages: 

direct estimation of free surface and phase/ 
fluid interfaces; 

effective simulation of multiple moving bodies 
and particulate matter within fluid flows; 

direct estimation of Lagrangian derivatives 
(absence of non-linear convective terms in the 
balance equation Left-Hand Sides); 

effective numerical simulations of fast transient 
phenomena; 

no meshing; 

simple and non-iterative algorithm. 

On the other hand, SPH models are still affected by a 
few shortcomings, if compared with mesh-based CFD 
tools: 

higher computational time consumption at the 
same resolution (mean feature); 

limited accuracy, if considered as a general 
purpose CFD technique. 


SPH models are effective in several, but peculiar, 
application fields. Some of them are here briefly 
recalled: floods; sloshing (e.g. oscillatory flows within 
tanks, Souto-Iglesias et al. 2006, Delorme et al. 2009, 
Amicarelli et al. 2013); gravitational surface waves (e.g. 
Manenti et al. 2008, Patel et al. 2009, Antuono et al. 
2011, Liu et al. 2013, Omidvar et al. 2013, Farahani et al. 
2014); hydraulic turbines (Marongiu et al. 2009, 
Neuhauser et al. 2013); fast landslides (Kumar et al. 


2013); hydroelectric plant devices; bed load transport 
and erosion phenomena (Manenti et al., 2009 and 2012, 
Manenti 2011, Guandalini-Agate, 2011); liquid jets 
(Koukouvinis et al. 2013); pollutant dispersion; 
astrophysics (e.g. Price and Monaghan 2007; Bauer 
and Spingel 2012); magnetohydrodynamics (e.g. Price 
2012); multi-phase and multi-fluid flows (e.g. Kajtar- 
Monaghan, 2012). 2D SPH tsunami modelling 
(Shallow Water approximation) has been already 
tested, considering the ISEC 2004 Benchmark2, which 
reproduces Okushiri tsunami (Vacondio et al. 2012). 

However, it seems there is a general lack of numerical 
model validations on complex tsunami events and 
coastal floods, especially in 3D. 

In this context, this study is focused on the validation 
of the recent SPH model (Di Monaco et al., 2011) on a 
reference 3D tsunami test case (Benchmark 1-2009, 
ISEC). After this introduction, we recall the main 
features of the numerical model (Sec.2). We then 
describe the reference tsunami test case and the 
numerical configuration (Sec. 3). Further we show and 
discuss the results (Sec.4) and we finally report the 
overall conclusions (Sec.5). 


The SPH Numerical Model 


The numerical model used for this validation is based 
on the semi-analytic approach for complex boundary 
treatment (Vila, 1999). Its basic features are deeply 
described in Di Monaco et al. (2011) and here briefly 
reported. 


Navier-Stokes equation and the continuity equation 
are first considered in the continuum: 


dUj 

dt 



d 2 u, 


p dp 

— T S a g+V 

p dx t dx 


J = 1,2,3 


( 3 ) 


dp 

dt 


= -pdivu 


where u is the velocity vector, p pressure, □ the fluid 
density, 8ij Kronecker's delta and t time (Einstein 
notation works for "j"). The pressure gradient term has 
been split after a simple product rule for derivatives, 
in order to be more easily treated. In fact these 
equations are then elaborated using the SPH particle 
approximation of (2) and computed at each fluid 
particle position, as described in the following. The 
viscosity term for compressible fluids can be neglected 
in Navier Stokes equation, even if we adopt a Weakly 
Compressible fluid approach. 


The SPH particle (discrete) approximation of a generic 
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derivative (2) can be computed at boundaries, 
according to the semi-analytic approach ("sa"): 


f-) =K/*-/o). 

dx i Isa b dx - 


dW„ 


dW„ 


^ + l(/-/o)^ 3 (4) 


6X: 


The summation is performed over all the fluid 
particles "b" (neighbouring particles of volume co) 
lying in the kernel support of the computational fluid 
particle ("o"). This formula also shows a boundary 
term, which is a convolution integral on the portion of 
the kernel support, lying outside the fluid domain. 
This portion then completes the kernel support, when 
it is truncated by boundaries. In this outer and 
fictitious volume (Vh'), we need to define a generic 
function f (pressure, velocity or density alternatively) 
to provide the proper boundary conditions. The semi- 
analytic approach, as developed by Di Monaco et al. 
(2011), follows this hypothesis: 

df I 


f = fsA + 


dx- 
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The functions and derivatives of these first order 
Taylor approximations are then assigned according to 
the following assumptions: 

' dp ' 


PSA = A)> 


PSA ~ A) * 


dp 

dx- 


(6) 


= 0 


SA 


These conditions try to represent a null normal 
gradient for the pressure displacement from hydrostatic 
conditions. Density is approximately considered to be 
a constant within the outer support Vh' (of the on- 
going computational particle). 


At the same time, free-slip conditions can be 
approximately set when estimating velocity in Vh. The 
velocity vector is assumed uniform in this region (still 
depending on the computational particle). It is 
decomposed in the sum of a vector normal to the 
boundary ( u SA n ) and the one tangential to it (u SA T ) . 

The first is represented by a linear extrapolation of the 
corresponding component of the computational fluid 
particle velocity. The latter is equal to the analogous 
value of the same particle (the subscript "w" refers to 
the local truncating frontier): 

USA = USAJ + —SA,n = 0 ~ A ) U 0.T + [(\ “ «0 ) ' 1w ] » w ( 7 ) 


In this case 4> s is zero (free-slip conditions; for no-slip 
conditions 4> s =l). 


The continuity equation ("Weakly-Compressible" SPH 
approach) can then be derived (Einstein's notation 
works for ""): 

dp o 


dt 


v ( 
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Finally the approximation of the momentum equation 
provides the following expression: 
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U r L dx i , 

vm represents the artificial viscosity (Monaghan, 2005), 
m the particle mass and r the relative distance between 
the neighbouring particle/boundary and the computa- 
tional particle. 


The second and third terms in the Right-Hand Side of 
(9) represent the inner and boundary terms related to 
pressure gradients. The next summation (Basa et al. 
2009) and integral represent the molecular viscous 
terms (the latter being proposed only for no-slip 
conditions by Di Monaco et al. 2011, but not yet 
validated in 3D). We finally have two artificial 
viscosity terms, commonly used to stabilize SPH 
models and defined by Monaghan (2005, inner term) 
and Di Monaco et al. 2011 (boundary term). They are 
here activated, only for computational particles 
approaching neighbouring particles or boundaries, 
according to Di Monaco et al. (2011). 


Time integration is performed with a staggered first 
order Euler scheme, where the continuity equation is 
translated in time of dt/2 with respect to the 
momentum equation (Di Monaco et al. 2011, Violeau- 
Leroy 2014). 


Finally a barotropic equation of state is linearized for 
small density oscillations as follows (cref and pref are the 
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reference sound speed and density): 

P = C )ef{P~ Pref ) ( 10 ) 

The 3D Tsunami Test Case 

The reference benchmark was proposed by the Johns 
Hopkins University and belongs to ISEC database. The 
3D domain and monitoring points are shown in FIG 1. 



FIG. 1 TOP: 3D DOMAIN, WAVEMAKER AND VELOCITY 
MONITORING POINTS (Vl-3). CENTRE: SEA LEVEL 
MONITORING LINES (Pl-12). BOTTOM: HEIGHT OF THE FLUID- 
BOTTOM INTERFACE (batimetry). THE PADDLE (wavemaker) IS 
SIMULATED USING SPH PARTICLES WITH IMPOSED 
PARAMETERS. THE VERTICAL AXIS IS STRETCHED BY A 
FACTOR OF 2, IN ALL THE FIGURES. 



A solitary wave (or tsunami) is generated on the left of 
the domain and propagates over a pyramidal shaped 
inclined shelf. In particular, a wave generator (left-side 
of the domain) produces a wave of amplitude h*=0.39 
m over the quiet free surface elevation (z=0.78 m; the 
origin of the vertical axis is provided by the lowest 
point in the domain). The tsunami propagates towards 
the coastline (right-side of the domain). The irregular 
shelf is represented by a flat bottom, which evolves 
into an inclined shelf close to the the coast, at whose 
centre a pyramidal structure rises up. The overall 
phenomenon lasts 45s, but the most dynamics refers to 
the first 15-20s. 

The numerical simulation has been configured 
accordingly to the experimental intial and boundary 
conditions. The parallelepiped inscribing the 
numerical domain is approximately 70x45x3 m 3 sized. 
The ratio between the kernel length size (h) and the 
particle size (dx) is equal to h/dx=1.25, with 
approximately dx=0.100 m. Time integration is 
performed choosing a high CFL (Courant-Friedrichs- 
Lewy) number of around 0.7, using a staggered Euler 
scheme. The wave generator movement is imposed, 
according to the experimental dataset. The lateral 
boundaries are fixed (semi-analytical approach) and 
water can move beyond the coast-line. 

Results 

The model results are here presented in terms of 
pressure and velocity fields, together with validations 
of the numerical sea level and velocity time series. 



Velocity (m/s) 

* * 02 0 .4 06 qa 1 1-2 14 u 

Time: 3.50 s q 



FIG. 2 3D FIELDS OF THE ABSOLUTE VALUE OF VELOCITY (top views). 
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FIG. 3 SEQUENCE OF PRESSURE 3D FIELDS (vertical slice through the symmetry plane of the promontory). 
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FIG. 4 SEQUENCE OF 3D FIELDS OF THE ABSOLUTE VALUE OF VELOCITY (vertical slice through the symmetry plane of the promontory). 


FIG 2 reports a sequence of velocity fields (top views). 
The wave generation lasts around 3s. The first stage of 
the wave propagation is mainly 2D (t=3.5s). When the 
tsunami reaches the submerged promontory (t=5.5s), 
the wave group speed has a local minimum (lowest 
bottom depth). It gradually increases towards the 
lateral boundaries of the domain, with an increasing 
bottom depth (at the same time the surface fluid 
velocity is highest over the promontory and decreseas 
far away). The wave direction locally converges to the 
promontory axis (wave refraction) while the wave 
begins to break (t around 6.5s). The phenomenon is 
affected by wave reflections both at the left and right 
domain boundaries. Wave interactions, minor 
stationary structures and jumps locally take place. 


until water reaches again an almost quiet equilibrium 
(after t=30-35s). 

The 3D fields of pressure and velocity are also 
described by a sequence of vertical slices passing for 
the promontory bary centres (FIG 3 and FIG 4). 

Pressure can locally fluctuate around the local 
hydrostatic vertical profiles (FIG 3). Despite some 
small-scale anomalies and a very few fluid particles 
penetrating the wavemaker, the model represents a 
regular pressure field and the limited errors do not 
propagate in time. Moreover we are not interested in 
very accurate estimations of pressure forces at 
boundaries, as this study focuses on the assessment of 
sea levels and velocity during a tsunami event. 
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The sequence of fields in FIG 4 shows the time and 
space distribution of velocity, whose maximum values 
reach around 1.6-1.8m/s. Velocity decreases from the 
free surface to the bottom, or moving far away from 
the solitary wave crests. High velocities are simulated 
in the breaking zone, then close to the paddle (during 
the wave generation) and beyond the coastline, during 
the flood. This event can be clearly detected at t=19.5s, 
with water covering more than a tenth of meters 
beyond the original coastline. Some details of the wave 
breaking are clearly represented in FIG. 5, which refers 
to a further simulation, with a slightly higher 
resolution (dx=0.08m) than the reference simulation. 



FIG. 5 Details of the wave breaking (dx=0.08m). 

Validations are performed against measures from 15 
monitoring points or lines. All the plots are 
normalized using reference height (h*=0.390m), 
velocity (u* = -J2gh* = 2.76m j and time scale t*(=h*/u* 

=0.141s), so that we treat non-dimensional variables 
(H=h/h*, U=u/u*, T=t/t*). 

In particular, velocity estimations are validated in 
terms of time series of the x-component (u) of the 
mean velocity, at three monitoring points (Vl-3). 

The monitor VI is placed along the vertical line 
crossing the vertex of the promontory (FIG 1). Its 
position, expressed in meters, is (13.00; 0.00; 0.75 m). 
The main forward and backward movements of water 
occur between T=32 and T=70. In this period we 
register the passage of the solitary wave, with water 
particles describing elliptical-like non-stationary 
trajectories. After this period a secondary wave is 
detected. The simulated values correctly reproduce the 
main features of u time evolution, even if some 
underestimations affect the amplitude of the peaks. 

The velocity monitors V2 and V3 (FIG 1, centre and 
bottom, respectively) are closer to the coastline than 
VI and transversally displaced with respect to the 
promontory. In these cases the measured time series 
show some lack of data (reference: ISEC web site). 
However several comparisons can still be performed. 
They show a satisfactory agreement between 
simulated and measured values. Nevertheless we can 
not catch, at this space resolution, the sudden sign 
inversion, which occurs at the monitor V2 (T=130). 




FIG. 6 TIME SERIES OF THE X-COMPONENT OF THE MEAN 
VELOCITY AT THE MONITORING POINTS Vl-3. COMPARISONS 
BETWEEN THE SIMULATED AND MEASURED VALUES. 
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FIG. 7 SEA LEVEL TIME SERIES AT THE MONITORING POINTS 
Pll, P9, P7, P5, P3. COMPARISONS BETWEEN THE SIMULATED 
AND MEASURED VALUES. 

Two asymmetric rows of measuring probes (Pl-12, 
FIG 1) are used to validate the estimation of the sea 
levels (FIG 7, FIG 8, FIG 9, FIG 10). All the curves 
clearly detect the passage of the solitary wave, which 
is followed by a secondary wave motion. 
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FIG. 8 SEA LEVEL TIME SERIES AT THE MONITORING POINTS 
P12, P10, P8, P6, P4. COMPARISONS BETWEEN THE SIMULATED 
AND MEASURED VALUES. 

The increment in sea level is accompanied by an 
increase in the wave group velocity and the 
fundamental frequency. The normalized sea levels 
show its maximum around the promontory (refraction 
effects), as already stated (surface wave velocity goes 
with the square root of the fluid depth, under shallow 
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water approximations, so that power fluxes and sea 
level peaks tend to concentrate around promontories). 
The comparisons between all the simulated values and 
the corresponding measures show the good reliability 
of the model in estimating sea levels. Nevertheless we 
notice some lack of numerical data (a numerical probe 
invalidates some sea level values when the particle 
distribution around it is not considered enough 
populated) and some peak underestimations, 
occurring at this space resolution. We can also notice 
that we slightly overestimate the sea level after the 
peak, especially at probe P3. This could involve an 
underestimation of the fluid velocity in the breaking 
zone or an overestimation of the wave power. 
However, the first hypothesis seems being invalidated 
by the velocity comparisons (FIG 6) which do not 
show any underestimation in the breaking zone. 
Finally, we underline the importance of a correct 
evaluation of the velocity field down-flow the wave 
breaking to assess the tsunami destructive power and 
run up on land. 



FIG. 9 SEA LEVEL TIME SERIES AT THE MONITORING POINT PI. 
COMPARISONS BETWEEN THE SIMULATED AND 
MEASURED VALUES. 



FIG. 10 SEA LEVEL TIME SERIES AT THE MONITORING POINT 
P2. COMPARISONS BETWEEN THE SIMULATED AND 
MEASURED VALUES. 

Conclusions 

A 3D SPH numerical model has been validated over a 


laboratory tsunami test case. The results show a good 
reliability of the model in representing the 
propagation of a solitary wave on a 3D shelf, both in 
terms of velocity and sea level. 

The application field of this kind of validations can 
involve gravitational surface wave propagation in 
natural and artificial basins, the corresponding 
interactions with complex 3D shelves and structures, 
together with the associated coastal flood events. 
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